Code
library(tidyverse) # Includes dplyr, ggplot2, readr, etc.
library(sjPlot) # to set and configure plotting theme
library(rnaturalearth) # to create base world map with countriesWe load the necessary libraries
library(tidyverse) # Includes dplyr, ggplot2, readr, etc.
library(sjPlot) # to set and configure plotting theme
library(rnaturalearth) # to create base world map with countriesWe configure the plotting theme based on the data visualization blog from Cédric Scherer: https://www.cedricscherer.com/
theme_Scherer <- function (base_size = 12, base_family = "Helvetica") {
half_line <- base_size/2
theme(
line = element_line(color = "black", linewidth = .5,
linetype = 1, lineend = "butt"),
rect = element_rect(fill = "white", color = "black",
linewidth = .5, linetype = 1),
text = element_text(family = base_family, face = "plain",
color = "black", size = base_size,
lineheight = .9, hjust = .5, vjust = .5,
angle = 0, margin = margin(), debug = FALSE),
axis.line = element_blank(),
axis.line.x = NULL,
axis.line.y = NULL,
axis.text = element_text(size = base_size * 1.1, color = "gray30"),
axis.text.x = element_text(margin = margin(t = .8 * half_line/2),
vjust = 1),
axis.text.x.top = element_text(margin = margin(b = .8 * half_line/2),
vjust = 0),
axis.text.y = element_text(margin = margin(r = .8 * half_line/2),
hjust = 1),
axis.text.y.right = element_text(margin = margin(l = .8 * half_line/2),
hjust = 0),
axis.ticks = element_line(color = "gray30", linewidth = .7),
axis.ticks.length = unit(half_line / 1.5, "pt"),
axis.ticks.length.x = NULL,
axis.ticks.length.x.top = NULL,
axis.ticks.length.x.bottom = NULL,
axis.ticks.length.y = NULL,
axis.ticks.length.y.left = NULL,
axis.ticks.length.y.right = NULL,
axis.title.x = element_text(margin = margin(t = half_line),
vjust = 1, size = base_size * 1.3,
face = "bold"),
axis.title.x.top = element_text(margin = margin(b = half_line),
vjust = 0),
axis.title.y = element_text(angle = 90, vjust = 1,
margin = margin(r = half_line),
size = base_size * 1.3, face = "bold"),
axis.title.y.right = element_text(angle = -90, vjust = 0,
margin = margin(l = half_line)),
legend.background = element_rect(color = NA),
legend.spacing = unit(.4, "cm"),
legend.spacing.x = NULL,
legend.spacing.y = NULL,
legend.margin = margin(.2, .2, .2, .2, "cm"),
legend.key = element_rect(fill = "gray95", color = "white"),
legend.key.size = unit(1.2, "lines"),
legend.key.height = NULL,
legend.key.width = NULL,
legend.text = element_text(size = rel(.8)),
legend.text.align = NULL,
legend.title = element_text(hjust = 0),
legend.title.align = NULL,
legend.position = "right",
legend.direction = NULL,
legend.justification = "center",
legend.box = NULL,
legend.box.margin = margin(0, 0, 0, 0, "cm"),
legend.box.background = element_blank(),
legend.box.spacing = unit(.4, "cm"),
panel.background = element_rect(fill = "white", color = NA),
panel.border = element_rect(color = "gray30",
fill = NA, linewidth = .7),
panel.grid.major = element_line(color = "gray90", linewidth = 1),
panel.grid.minor = element_line(color = "gray90", linewidth = .5,
linetype = "dashed"),
panel.spacing = unit(base_size, "pt"),
panel.spacing.x = NULL,
panel.spacing.y = NULL,
panel.ontop = FALSE,
strip.background = element_rect(fill = "white", color = "gray30"),
strip.text = element_text(color = "black", size = base_size),
strip.text.x = element_text(margin = margin(t = half_line,
b = half_line)),
strip.text.y = element_text(angle = -90,
margin = margin(l = half_line,
r = half_line)),
strip.text.y.left = element_text(angle = 90),
strip.placement = "inside",
strip.placement.x = NULL,
strip.placement.y = NULL,
strip.switch.pad.grid = unit(0.1, "cm"),
strip.switch.pad.wrap = unit(0.1, "cm"),
plot.background = element_rect(color = NA),
plot.title = element_text(size = base_size * 1.8, hjust = .5,
vjust = 1, face = "bold",
margin = margin(b = half_line * 1.2)),
plot.title.position = "panel",
plot.subtitle = element_text(size = base_size * 1.3,
hjust = .5, vjust = 1,
margin = margin(b = half_line * .9)),
plot.caption = element_text(size = rel(0.9), hjust = 1, vjust = 1,
margin = margin(t = half_line * .9)),
plot.caption.position = "panel",
plot.tag = element_text(size = rel(1.2), hjust = .5, vjust = .5),
plot.tag.position = "topleft",
plot.margin = margin(rep(base_size, 4)),
complete = TRUE
)
}
set_theme(theme_Scherer())Then we can load our data:
base_path <- "C:/Users/kobed/Documents/TOAD-PIKlegacy" # Where is the code repo stored?
# List of crops
crops <- c("mai", "wwh", "ri1", "soy", "aggr")
# Read each .rds file into a named list
crop_data <- lapply(crops, function(crop) {
readRDS(file.path(base_path, paste0("extremes_", crop, ".rds")))
})
names(crop_data) <- cropsFirst we prepare the data
dat_long <- crop_data[["aggr"]] %>%
dplyr::select(lat, lon, year, divtrend_sim, divtrend_obs, model, climate_extreme) %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
divtrend_obs_perc = (divtrend_obs - 1)*100
) %>%
pivot_longer(
cols = c(divtrend_obs_perc, divtrend_sim_perc), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "divtrend"
) %>%
mutate(
model = case_when(
type == "divtrend_obs_perc" ~ "benchmark",
type == "divtrend_sim_perc" ~ model
)
) %>%
distinct() Then we can calculate the proportion of underestimation and the median detrended yield for the benchmark data
underestimation_rate <- crop_data[["aggr"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
benchmark_medians <- dat_long %>%
filter(model == "benchmark") %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(divtrend, na.rm = TRUE))Prepare the plot by selecting the name for the benchmark and creating a fixed order to list the models
# Define the name of the benchmark distribution
benchmark_name <- "benchmark"
# Model order
order <- rev(c(
"benchmark",
"ensemble",
"acea",
"simplace-lintul5",
"dssat-pythia",
"pepic",
"cygma1p74",
"ldndc",
"pdssat",
"epic-iiasa",
"lpj-guess",
"promet",
"lpjml",
"isam",
"crover"
))Now we can create the raincloud plot
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
filter(is.finite(divtrend), abs(divtrend) < 1000) %>% # Filter out Inf values and keep within range
ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add black vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
x = "Yield anomalies in %",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "fig4_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)dat_long <- crop_data[["ri1"]] %>%
dplyr::select(lat, lon, year, divtrend_sim, divtrend_obs, model, climate_extreme) %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
divtrend_obs_perc = (divtrend_obs - 1)*100
) %>%
pivot_longer(
cols = c(divtrend_obs_perc, divtrend_sim_perc), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "divtrend"
) %>%
mutate(
model = case_when(
type == "divtrend_obs_perc" ~ "benchmark",
type == "divtrend_sim_perc" ~ model
)
) %>%
distinct()
underestimation_rate <- crop_data[["ri1"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
benchmark_medians <- dat_long %>%
filter(model == "benchmark") %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(divtrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
filter(is.finite(divtrend), abs(divtrend) < 1000) %>% # Filter out Inf values and keep within range
ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add black vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
x = "Yield anomalies in %",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig10_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)dat_long <- crop_data[["wwh"]] %>%
dplyr::select(lat, lon, year, divtrend_sim, divtrend_obs, model, climate_extreme) %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
divtrend_obs_perc = (divtrend_obs - 1)*100
) %>%
pivot_longer(
cols = c(divtrend_obs_perc, divtrend_sim_perc), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "divtrend"
) %>%
mutate(
model = case_when(
type == "divtrend_obs_perc" ~ "benchmark",
type == "divtrend_sim_perc" ~ model
)
) %>%
distinct()
underestimation_rate <- crop_data[["wwh"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
benchmark_medians <- dat_long %>%
filter(model == "benchmark") %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(divtrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
filter(is.finite(divtrend), abs(divtrend) < 1000) %>% # Filter out Inf values and keep within range
ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add black vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
x = "Yield anomalies in %",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig11_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)dat_long <- crop_data[["soy"]] %>%
dplyr::select(lat, lon, year, divtrend_sim, divtrend_obs, model, climate_extreme) %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
divtrend_obs_perc = (divtrend_obs - 1)*100
) %>%
pivot_longer(
cols = c(divtrend_obs_perc, divtrend_sim_perc), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "divtrend"
) %>%
mutate(
model = case_when(
type == "divtrend_obs_perc" ~ "benchmark",
type == "divtrend_sim_perc" ~ model
)
) %>%
distinct()
underestimation_rate <- crop_data[["soy"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
benchmark_medians <- dat_long %>%
filter(model == "benchmark") %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(divtrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
filter(is.finite(divtrend), abs(divtrend) < 1000) %>% # Filter out Inf values and keep within range
ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add black vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
x = "Yield anomalies in %",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig12_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)dat_long <- crop_data[["mai"]] %>%
dplyr::select(lat, lon, year, divtrend_sim, divtrend_obs, model, climate_extreme) %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
divtrend_obs_perc = (divtrend_obs - 1)*100
) %>%
pivot_longer(
cols = c(divtrend_obs_perc, divtrend_sim_perc), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "divtrend"
) %>%
mutate(
model = case_when(
type == "divtrend_obs_perc" ~ "benchmark",
type == "divtrend_sim_perc" ~ model
)
) %>%
distinct()
underestimation_rate <- crop_data[["mai"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
benchmark_medians <- dat_long %>%
filter(model == "benchmark") %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(divtrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
filter(is.finite(divtrend), abs(divtrend) < 1000) %>% # Filter out Inf values and keep within range
ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add black vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
x = "Yield anomalies in %",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig13_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)We need the unfiltered benchmark data to find out who the biggest producers are for this dataset.
# Read each .rds file into a named list
crop_datageneral <- lapply(crops, function(crop) {
readRDS(file.path(base_path, paste0("general_", crop, ".rds")))
})
general_sim_aggr <- readRDS(file.path(base_path, "general_sim_aggr.rds"))
names(crop_datageneral) <- cropsFirst we calculate who the 5 largest producing countries are and also calculate how much tons they produce and what the global share is (percentage)
crop_datageneral[["aggr"]] %>%
left_join(general_sim_aggr, by = c("lon", "lat", "year")) %>%
mutate(prod = yield * total_area) %>%
group_by(ctr) %>%
summarise(ctr_prod = sum(prod, na.rm = TRUE), .groups = "drop") %>%
mutate(global_prod = sum(ctr_prod),
global_share = (ctr_prod / global_prod) *100) %>%
slice_max(order_by = ctr_prod, n = 5) %>%
dplyr::select(ctr, ctr_prod, global_share)# A tibble: 5 × 3
ctr ctr_prod global_share
<chr> <dbl> <dbl>
1 CHN 5.63e15 23.5
2 USA 5.53e15 23.1
3 BRA 1.71e15 7.14
4 ARG 1.05e15 4.39
5 IND 7.15e14 2.99
Then we plot
underestimation_rate <- crop_data[["aggr"]] %>%
filter(ctr == "BRA" |ctr == "USA" | ctr == "IND" | ctr == "ARG" | ctr == "CHN") %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
dat_long <- crop_data[["aggr"]] %>%
filter(ctr == "BRA" |ctr == "USA" | ctr == "IND" | ctr == "ARG" | ctr == "CHN") %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
divtrend_obs_perc = (divtrend_obs - 1)*100
) %>%
dplyr::select(lat, lon, year, divtrend_sim_perc, divtrend_obs_perc, model, climate_extreme) %>%
pivot_longer(
cols = c(divtrend_obs_perc, divtrend_sim_perc), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "divtrend"
) %>%
mutate(
model = case_when(
type == "divtrend_obs_perc" ~ "benchmark",
type == "divtrend_sim_perc" ~ model
)
) %>%
distinct()
# Calculate the median yield anomaly for the benchmark model for each climate extreme
benchmark_medians <- dat_long %>%
filter(model == benchmark_name) %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(divtrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
filter(is.finite(divtrend), abs(divtrend) < 1000) %>% # Filter out Inf values and keep within range
ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 2, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add red vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
subtitle = "Top 5 producers: CHN (23.5%), USA (23.1%), BRA (7.14%), ARG(4.39%), IND(2.99%)",
x = "Yield anomalies in %",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig14_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)First we calculate who the 5 largest producing countries are and also calculate how much tons they produce and what the global share is (percentage)
crop_datageneral[["ri1"]] %>%
mutate(prod = yield * total_area) %>%
group_by(ctr) %>%
summarise(ctr_prod = sum(prod, na.rm = TRUE), .groups = "drop") %>%
mutate(global_prod = sum(ctr_prod),
global_share = (ctr_prod / global_prod) *100) %>%
slice_max(order_by = ctr_prod, n = 5) %>%
dplyr::select(ctr, ctr_prod, global_share)# A tibble: 5 × 3
ctr ctr_prod global_share
<chr> <dbl> <dbl>
1 CHN 8.90e14 30.7
2 IND 4.11e14 14.2
3 IDN 2.02e14 6.98
4 THA 1.95e14 6.71
5 BRA 1.84e14 6.35
Then we plot
underestimation_rate <- crop_data[["ri1"]] %>%
filter(ctr == "BRA" |ctr == "THA" | ctr == "IND" | ctr == "IDN" | ctr == "CHN") %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
dat_long <- crop_data[["ri1"]] %>%
filter(ctr == "BRA" |ctr == "THA" | ctr == "IND" | ctr == "IDN" | ctr == "CHN") %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
divtrend_obs_perc = (divtrend_obs - 1)*100
) %>%
dplyr::select(lat, lon, year, divtrend_sim_perc, divtrend_obs_perc, model, climate_extreme) %>%
pivot_longer(
cols = c(divtrend_obs_perc, divtrend_sim_perc), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "divtrend"
) %>%
mutate(
model = case_when(
type == "divtrend_obs_perc" ~ "benchmark",
type == "divtrend_sim_perc" ~ model
)
) %>%
distinct()
# Calculate the median yield anomaly for the benchmark model for each climate extreme
benchmark_medians <- dat_long %>%
filter(model == benchmark_name) %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(divtrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
filter(is.finite(divtrend), abs(divtrend) < 1000) %>% # Filter out Inf values and keep within range
ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 2, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add red vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
subtitle = "Top 5 producers: CHN (30.7%), IND (14.2%), IDN (6.98%), THA(6.71%), BRA(6.35%)",
x = "Yield anomalies in %",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig15_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)First we calculate who the 5 largest producing countries are and also calculate how much tons they produce and what the global share is (percentage)
crop_datageneral[["wwh"]] %>%
mutate(prod = yield * total_area) %>%
group_by(ctr) %>%
summarise(ctr_prod = sum(prod, na.rm = TRUE), .groups = "drop") %>%
mutate(global_prod = sum(ctr_prod),
global_share = (ctr_prod / global_prod) *100) %>%
slice_max(order_by = ctr_prod, n = 5) %>%
dplyr::select(ctr, ctr_prod, global_share)# A tibble: 5 × 3
ctr ctr_prod global_share
<chr> <dbl> <dbl>
1 USA 4.15e14 14.1
2 FRA 3.75e14 12.7
3 CHN 3.35e14 11.4
4 DEU 2.51e14 8.53
5 TUR 1.76e14 5.97
Then we plot
underestimation_rate <- crop_data[["wwh"]] %>%
filter(ctr == "FRA" |ctr == "USA" | ctr == "DEU" | ctr == "TUR" | ctr == "CHN") %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
dat_long <- crop_data[["wwh"]] %>%
filter(ctr == "FRA" |ctr == "USA" | ctr == "DEU" | ctr == "TUR" | ctr == "CHN") %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
divtrend_obs_perc = (divtrend_obs - 1)*100
) %>%
dplyr::select(lat, lon, year, divtrend_sim_perc, divtrend_obs_perc, model, climate_extreme) %>%
pivot_longer(
cols = c(divtrend_obs_perc, divtrend_sim_perc), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "divtrend"
) %>%
mutate(
model = case_when(
type == "divtrend_obs_perc" ~ "benchmark",
type == "divtrend_sim_perc" ~ model
)
) %>%
distinct()
# Calculate the median yield anomaly for the benchmark model for each climate extreme
benchmark_medians <- dat_long %>%
filter(model == benchmark_name) %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(divtrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
filter(is.finite(divtrend), abs(divtrend) < 1000) %>% # Filter out Inf values and keep within range
ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 2, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add red vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
subtitle = "Top 5 producers: USA (14.1%), FRA (12.7%), CHN (11.4%), DEU(8.53%), TUR(5.97%)",
x = "Yield anomalies in %",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig16_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)First we calculate who the 5 largest producing countries are and also calculate how much tons they produce and what the global share is (percentage)
crop_datageneral[["soy"]] %>%
mutate(prod = yield * total_area) %>%
group_by(ctr) %>%
summarise(ctr_prod = sum(prod, na.rm = TRUE), .groups = "drop") %>%
mutate(global_prod = sum(ctr_prod),
global_share = (ctr_prod / global_prod) *100) %>%
slice_max(order_by = ctr_prod, n = 5) %>%
dplyr::select(ctr, ctr_prod, global_share)# A tibble: 5 × 3
ctr ctr_prod global_share
<chr> <dbl> <dbl>
1 USA 5.92e14 41.4
2 BRA 4.04e14 28.2
3 ARG 2.40e14 16.8
4 CHN 1.11e14 7.78
5 IND 1.96e13 1.37
Then we plot
underestimation_rate <- crop_data[["soy"]] %>%
filter(ctr == "BRA" |ctr == "USA" | ctr == "IND" | ctr == "ARG" | ctr == "CHN") %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
dat_long <- crop_data[["soy"]] %>%
filter(ctr == "BRA" |ctr == "USA" | ctr == "IND" | ctr == "ARG" | ctr == "CHN") %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
divtrend_obs_perc = (divtrend_obs - 1)*100
) %>%
dplyr::select(lat, lon, year, divtrend_sim_perc, divtrend_obs_perc, model, climate_extreme) %>%
pivot_longer(
cols = c(divtrend_obs_perc, divtrend_sim_perc), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "divtrend"
) %>%
mutate(
model = case_when(
type == "divtrend_obs_perc" ~ "benchmark",
type == "divtrend_sim_perc" ~ model
)
) %>%
distinct()
# Calculate the median yield anomaly for the benchmark model for each climate extreme
benchmark_medians <- dat_long %>%
filter(model == benchmark_name) %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(divtrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
filter(is.finite(divtrend), abs(divtrend) < 1000) %>% # Filter out Inf values and keep within range
ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 2, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add red vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
subtitle = "Top 5 producers: USA (41.4%), BRA (28.2%), ARG (16.8%), CHN(7.78%), IND(1.37%)",
x = "Yield anomalies in %",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig17_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)First we calculate who the 5 largest producing countries are and also calculate how much tons they produce and what the global share is (percentage)
crop_datageneral[["mai"]] %>%
mutate(prod = yield * total_area) %>%
group_by(ctr) %>%
summarise(ctr_prod = sum(prod, na.rm = TRUE), .groups = "drop") %>%
mutate(global_prod = sum(ctr_prod),
global_share = (ctr_prod / global_prod) *100) %>%
slice_max(order_by = ctr_prod, n = 5) %>%
dplyr::select(ctr, ctr_prod, global_share)# A tibble: 5 × 3
ctr ctr_prod global_share
<chr> <dbl> <dbl>
1 USA 2.30e15 40.6
2 CHN 9.30e14 16.4
3 BRA 4.49e14 7.91
4 ARG 2.47e14 4.35
5 FRA 2.41e14 4.24
Then we plot
underestimation_rate <- crop_data[["mai"]] %>%
filter(ctr == "BRA" |ctr == "USA" | ctr == "FRA" | ctr == "ARG" | ctr == "CHN") %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
dat_long <- crop_data[["mai"]] %>%
filter(ctr == "BRA" |ctr == "USA" | ctr == "FRA" | ctr == "ARG" | ctr == "CHN") %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
divtrend_obs_perc = (divtrend_obs - 1)*100
) %>%
dplyr::select(lat, lon, year, divtrend_sim_perc, divtrend_obs_perc, model, climate_extreme) %>%
pivot_longer(
cols = c(divtrend_obs_perc, divtrend_sim_perc), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "divtrend"
) %>%
mutate(
model = case_when(
type == "divtrend_obs_perc" ~ "benchmark",
type == "divtrend_sim_perc" ~ model
)
) %>%
distinct()
# Calculate the median yield anomaly for the benchmark model for each climate extreme
benchmark_medians <- dat_long %>%
filter(model == benchmark_name) %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(divtrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
filter(is.finite(divtrend), abs(divtrend) < 1000) %>% # Filter out Inf values and keep within range
ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 2, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add red vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
subtitle = "Top 5 producers: USA (40.6%), CHN (16.4%), BRA (7.91%), ARG(4.35%), FRA(4.24%)",
x = "Yield anomalies in %",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig18_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)dat_long <- crop_data[["aggr"]] %>%
dplyr::select(lat, lon, year, difftrend_sim, difftrend_obs, model, climate_extreme) %>%
pivot_longer(
cols = c(difftrend_obs, difftrend_sim), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "difftrend"
) %>%
mutate(
model = case_when(
type == "difftrend_obs" ~ "benchmark",
type == "difftrend_sim" ~ model
)
) %>%
distinct()
underestimation_rate <- crop_data[["aggr"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error < 0) * 100) # Calculate % of values below 0
benchmark_medians <- dat_long %>%
filter(model == "benchmark") %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(difftrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = difftrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add black vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
x = "Yield anomalies in tons/ha",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-3, 3)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig19_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)dat_long <- crop_data[["ri1"]] %>%
dplyr::select(lat, lon, year, difftrend_sim, difftrend_obs, model, climate_extreme) %>%
pivot_longer(
cols = c(difftrend_obs, difftrend_sim), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "difftrend"
) %>%
mutate(
model = case_when(
type == "difftrend_obs" ~ "benchmark",
type == "difftrend_sim" ~ model
)
) %>%
distinct()
underestimation_rate <- crop_data[["ri1"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error < 0) * 100) # Calculate % of values below 0
benchmark_medians <- dat_long %>%
filter(model == "benchmark") %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(difftrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = difftrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add black vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
x = "Yield anomalies in tons/ha",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-3, 3)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig20_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)dat_long <- crop_data[["wwh"]] %>%
dplyr::select(lat, lon, year, difftrend_sim, difftrend_obs, model, climate_extreme) %>%
pivot_longer(
cols = c(difftrend_obs, difftrend_sim), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "difftrend"
) %>%
mutate(
model = case_when(
type == "difftrend_obs" ~ "benchmark",
type == "difftrend_sim" ~ model
)
) %>%
distinct()
underestimation_rate <- crop_data[["wwh"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error < 0) * 100) # Calculate % of values below 0
benchmark_medians <- dat_long %>%
filter(model == "benchmark") %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(difftrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = difftrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add black vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
x = "Yield anomalies in tons/ha",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-3, 3)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig21_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)dat_long <- crop_data[["soy"]] %>%
dplyr::select(lat, lon, year, difftrend_sim, difftrend_obs, model, climate_extreme) %>%
pivot_longer(
cols = c(difftrend_obs, difftrend_sim), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "difftrend"
) %>%
mutate(
model = case_when(
type == "difftrend_obs" ~ "benchmark",
type == "difftrend_sim" ~ model
)
) %>%
distinct()
underestimation_rate <- crop_data[["soy"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error < 0) * 100) # Calculate % of values below 0
benchmark_medians <- dat_long %>%
filter(model == "benchmark") %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(difftrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = difftrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add black vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
x = "Yield anomalies in tons/ha",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-3, 3)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig22_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)dat_long <- crop_data[["mai"]] %>%
dplyr::select(lat, lon, year, difftrend_sim, difftrend_obs, model, climate_extreme) %>%
pivot_longer(
cols = c(difftrend_obs, difftrend_sim), # Adjust to the actual column names you want to pivot
names_to = "type",
values_to = "difftrend"
) %>%
mutate(
model = case_when(
type == "difftrend_obs" ~ "benchmark",
type == "difftrend_sim" ~ model
)
) %>%
distinct()
underestimation_rate <- crop_data[["mai"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error < 0) * 100) # Calculate % of values below 0
benchmark_medians <- dat_long %>%
filter(model == "benchmark") %>%
group_by(climate_extreme) %>%
summarise(median_yield = median(difftrend, na.rm = TRUE))
# Create the raincloud plot
raincloud_perc <- dat_long %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = difftrend, y = model, fill = model == benchmark_name)) +
# Add half-eye plot
ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
scale_fill_manual(values = c("TRUE" = "#00A087FF",
"FALSE" = "#7E6148FF"))+
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add black vertical line for the benchmark median
geom_vline(data = benchmark_medians, aes(xintercept = median_yield),
color = "black", size = 2) +
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, color = "black") + # Adjust text size and positioning
# Labels and annotations
labs(title = "",
x = "Yield anomalies in tons/ha",
y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-3, 3)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_perc# Save the figure
ggsave(file.path(base_path, "suppfig23_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)Instead of separating benchmark from simulation data we can also create the raincloud plots for the model errors (benchmark - simulation data points)
underestimation_rate <- crop_data[["aggr"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["aggr"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
filter(is.finite(error_perc), abs(error_perc) < 1000) %>% # Filter out Inf values and keep within range
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = error_perc, y = model)) +
coord_trans(x = scales::modulus_trans(p = 0.5)) + # Adjust p for stronger compression
# Add half-eye plot with fill based on whether the error is negative or positive
ggdist::stat_halfeye(aes(fill = stat(x < 0)), # Logical fill based on error sign
adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add a vertical line at 0 (benchmark median)
geom_vline(aes(xintercept = 0), color = "black", size = 2) +
# Custom color scale where negative values are red and positive values are blue
scale_fill_manual(values = c("grey", "#DC0000FF") )+
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = -50.25, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") + # Adjust text size and positioning
# Labels and annotation
labs(x = "Model Error (benchmark - simulated) in % yield anomalies", y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_percggsave(file.path(base_path, "suppfig24_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)underestimation_rate <- crop_data[["ri1"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["ri1"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
filter(is.finite(error_perc), abs(error_perc) < 1000) %>% # Filter out Inf values and keep within range
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = error_perc, y = model)) +
coord_trans(x = scales::modulus_trans(p = 0.5)) + # Adjust p for stronger compression
# Add half-eye plot with fill based on whether the error is negative or positive
ggdist::stat_halfeye(aes(fill = stat(x < 0)), # Logical fill based on error sign
adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add a vertical line at 0 (benchmark median)
geom_vline(aes(xintercept = 0), color = "black", size = 2) +
# Custom color scale where negative values are red and positive values are blue
scale_fill_manual(values = c("grey", "#DC0000FF") )+
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = -50.25, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") + # Adjust text size and positioning
# Labels and annotation
labs(x = "Model Error (benchmark - simulated) in % yield anomalies", y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_percggsave(file.path(base_path, "suppfig25_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)underestimation_rate <- crop_data[["wwh"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["wwh"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
filter(is.finite(error_perc), abs(error_perc) < 1000) %>% # Filter out Inf values and keep within range
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = error_perc, y = model)) +
coord_trans(x = scales::modulus_trans(p = 0.5)) + # Adjust p for stronger compression
# Add half-eye plot with fill based on whether the error is negative or positive
ggdist::stat_halfeye(aes(fill = stat(x < 0)), # Logical fill based on error sign
adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add a vertical line at 0 (benchmark median)
geom_vline(aes(xintercept = 0), color = "black", size = 2) +
# Custom color scale where negative values are red and positive values are blue
scale_fill_manual(values = c("grey", "#DC0000FF") )+
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = -50.25, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") + # Adjust text size and positioning
# Labels and annotation
labs(x = "Model Error (benchmark - simulated) in % yield anomalies", y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_percggsave(file.path(base_path, "suppfig26_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)underestimation_rate <- crop_data[["soy"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["soy"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
filter(is.finite(error_perc), abs(error_perc) < 1000) %>% # Filter out Inf values and keep within range
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = error_perc, y = model)) +
coord_trans(x = scales::modulus_trans(p = 0.5)) + # Adjust p for stronger compression
# Add half-eye plot with fill based on whether the error is negative or positive
ggdist::stat_halfeye(aes(fill = stat(x < 0)), # Logical fill based on error sign
adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add a vertical line at 0 (benchmark median)
geom_vline(aes(xintercept = 0), color = "black", size = 2) +
# Custom color scale where negative values are red and positive values are blue
scale_fill_manual(values = c("grey", "#DC0000FF") )+
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = -50.25, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") + # Adjust text size and positioning
# Labels and annotation
labs(x = "Model Error (benchmark - simulated) in % yield anomalies", y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_percggsave(file.path(base_path, "suppfig27_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)underestimation_rate <- crop_data[["mai"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error_perc < 0) * 100) # Calculate % of values below 0
# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["mai"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
filter(is.finite(error_perc), abs(error_perc) < 1000) %>% # Filter out Inf values and keep within range
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = error_perc, y = model)) +
coord_trans(x = scales::modulus_trans(p = 0.5)) + # Adjust p for stronger compression
# Add half-eye plot with fill based on whether the error is negative or positive
ggdist::stat_halfeye(aes(fill = stat(x < 0)), # Logical fill based on error sign
adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add a vertical line at 0 (benchmark median)
geom_vline(aes(xintercept = 0), color = "black", size = 2) +
# Custom color scale where negative values are red and positive values are blue
scale_fill_manual(values = c("grey", "#DC0000FF") )+
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = -50.25, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") + # Adjust text size and positioning
# Labels and annotation
labs(x = "Model Error (benchmark - simulated) in % yield anomalies", y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-100, 100)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_percggsave(file.path(base_path, "suppfig28_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)underestimation_rate <- crop_data[["aggr"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error < 0) * 100) # Calculate % of values below 0
# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["aggr"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = error, y = model)) +
coord_trans(x = scales::modulus_trans(p = 0.5)) + # Adjust p for stronger compression
# Add half-eye plot with fill based on whether the error is negative or positive
ggdist::stat_halfeye(aes(fill = stat(x < 0)), # Logical fill based on error sign
adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add a vertical line at 0 (benchmark median)
geom_vline(aes(xintercept = 0), color = "black", size = 2) +
# Custom color scale where negative values are red and positive values are blue
scale_fill_manual(values = c("grey", "#DC0000FF") )+
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = -1.25, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") + # Adjust text size and positioning
# Labels and annotation
labs(x = "Model Error (benchmark - simulated) in tons/ha yield anomalies", y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-3, 3)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_percggsave(file.path(base_path, "suppfig29_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)underestimation_rate <- crop_data[["ri1"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error < 0) * 100) # Calculate % of values below 0
# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["ri1"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = error, y = model)) +
coord_trans(x = scales::modulus_trans(p = 0.5)) + # Adjust p for stronger compression
# Add half-eye plot with fill based on whether the error is negative or positive
ggdist::stat_halfeye(aes(fill = stat(x < 0)), # Logical fill based on error sign
adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add a vertical line at 0 (benchmark median)
geom_vline(aes(xintercept = 0), color = "black", size = 2) +
# Custom color scale where negative values are red and positive values are blue
scale_fill_manual(values = c("grey", "#DC0000FF") )+
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = -1.25, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") + # Adjust text size and positioning
# Labels and annotation
labs(x = "Model Error (benchmark - simulated) in tons/ha yield anomalies", y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-3, 3)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_percggsave(file.path(base_path, "suppfig30_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)underestimation_rate <- crop_data[["wwh"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error < 0) * 100) # Calculate % of values below 0
# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["wwh"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = error, y = model)) +
coord_trans(x = scales::modulus_trans(p = 0.5)) + # Adjust p for stronger compression
# Add half-eye plot with fill based on whether the error is negative or positive
ggdist::stat_halfeye(aes(fill = stat(x < 0)), # Logical fill based on error sign
adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add a vertical line at 0 (benchmark median)
geom_vline(aes(xintercept = 0), color = "black", size = 2) +
# Custom color scale where negative values are red and positive values are blue
scale_fill_manual(values = c("grey", "#DC0000FF") )+
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = -1.25, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") + # Adjust text size and positioning
# Labels and annotation
labs(x = "Model Error (benchmark - simulated) in tons/ha yield anomalies", y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-3, 3)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_percggsave(file.path(base_path, "suppfig31_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)underestimation_rate <- crop_data[["soy"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error < 0) * 100) # Calculate % of values below 0
# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["soy"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = error, y = model)) +
coord_trans(x = scales::modulus_trans(p = 0.5)) + # Adjust p for stronger compression
# Add half-eye plot with fill based on whether the error is negative or positive
ggdist::stat_halfeye(aes(fill = stat(x < 0)), # Logical fill based on error sign
adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add a vertical line at 0 (benchmark median)
geom_vline(aes(xintercept = 0), color = "black", size = 2) +
# Custom color scale where negative values are red and positive values are blue
scale_fill_manual(values = c("grey", "#DC0000FF") )+
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = -1.25, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") + # Adjust text size and positioning
# Labels and annotation
labs(x = "Model Error (benchmark - simulated) in tons/ha yield anomalies", y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-3, 3)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_percggsave(file.path(base_path, "suppfig32_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)underestimation_rate <- crop_data[["mai"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>% # Replace Inf with -1000
group_by(model, climate_extreme) %>%
summarise(under_rate = mean(error < 0) * 100) # Calculate % of values below 0
# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["mai"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
mutate(model = factor(model, levels = order)) %>% # Order the models using the defined order
ggplot(aes(x = error, y = model)) +
coord_trans(x = scales::modulus_trans(p = 0.5)) + # Adjust p for stronger compression
# Add half-eye plot with fill based on whether the error is negative or positive
ggdist::stat_halfeye(aes(fill = stat(x < 0)), # Logical fill based on error sign
adjust = 0.7,
width = 0.7, # Increase distribution width
justification = -0.1, # Move half-eye plot further from boxplot
.width = c(0.3, 0),
point_colour = "black",
alpha = 0.7) +
# Add boxplots (smaller width to give more space to distribution)
geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7,
fill = "white", color = "black", size = 1) +
# Add a vertical line at 0 (benchmark median)
geom_vline(aes(xintercept = 0), color = "black", size = 2) +
# Custom color scale where negative values are red and positive values are blue
scale_fill_manual(values = c("grey", "#DC0000FF") )+
# Add text annotations for the underestimation rate (place them to the right of the plot)
geom_text(data = underestimation_rate,
aes(x = -1.25, y = model, label = paste0(round(under_rate, 1), "%")), # Round percentage to 1 decimal
vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") + # Adjust text size and positioning
# Labels and annotation
labs(x = "Model Error (benchmark - simulated) in tons/ha yield anomalies", y = "") +
# Add facets for different climate extremes
facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
# Custom x-axis breaks
scale_x_continuous(breaks = seq(-5, 5, by = 1)) +
# Truncate the x-axis to focus on the range
coord_cartesian(xlim = c(-5, 5)) +
theme(legend.position = "none",
axis.title.x = element_text(size = 18), # Increase x-axis title size
axis.text.x = element_text(size = 18), # Increase x-axis text size
axis.text.y = element_text(size = 18), # Increase y-axis text size
strip.text = element_text(size = 18)) # Increase facet label size
raincloud_percggsave(file.path(base_path, "suppfig33_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)